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Abstract Wet granular materials in a quasi-static 
steady state shear flow have been studied with discrete 
particle simulations. Macroscopic quantities, consistent 
with the conservation laws of continuum theory, are ob¬ 
tained by time averaging and spatial coarse-graining. 
Initial studies involve understanding the effect of liq¬ 
uid content and liquid properties like the surface ten¬ 
sion on the macroscopic quantities. Two parameters 
of the liquid bridge contact model have been studied 
as the constitutive parameters that define the struc¬ 
ture of this model (i) the rupture distance of the liquid 
bridge model, which is proportional to the liquid con¬ 
tent, and (ii) the maximum adhesive force, as controlled 
by the surface tension of the liquid. Subsequently a 
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correlation is developed between these micro param¬ 
eters and the steady state cohesion in the limit of zero 
confining pressure. Furthermore, as second result, the 
macroscopic torque measured at the walls, which is an 
experimentally accessible parameter, is predicted from 
our simulation results as a dependence on the micro¬ 
parameters. Finally, the steady state cohesion of a real¬ 
istic non-linear liquid bridge contact model scales well 
with the steady state cohesion for a simpler linearized 
irreversible contact model with the same maximum ad¬ 
hesive force and equal energy dissipated per contact. 


1 Introduction 

Granular media are collections of microscopic grains 
having athermal interactions through dissipative, fric¬ 
tional or cohesive contact forces. External force leads to 
granular flow under the condition of applied shear stress 
exceeding the yield shear stress. After a finite shear 
strain, at constant rate, a steady state establishes with 
a typically lower shear stress, depending on both strain 
rate and pressure [1]. Most studies in granular physics 
focus on dry granular materials and their flow rheology. 
However, wet granular materials are ubiquitous in geol¬ 
ogy and many real world applications where interstitial 
liquid is present between the grains. Simplified models 
for capillary clusters [213] and wet granular gases jl] 
were introduced before. The rheology of flow for dense 
suspension of non-Brownian particles have been studied 
in Ref. [5l|6l[7]. We study the local rheology of weakly 
wetted granular materials in the quasistatic regime with 
the Discrete Element Method (DEM) using the open- 
source package MercuryDPM 0191 in a shear cell set¬ 
up, where the relative motion is confined to particles 
in a narrow region away from the walls, called shear 
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band unun]. We study partially saturated systems in 
the pendular regime, with a very low level of water con¬ 
tent, where the formation of liquid bridges between par¬ 
ticle pairs leads to development of microscopic tensile 
forces. This tensile forces generated at particle level re¬ 
sults in cohesion at macroscopic scale. Earlier studies 
have been done for liquid bridge in the pendular regime 
to understand the effect of liquid bridge volume and 
contact angle on different macroscopic quantities like 
the steady state cohesion, torque and shear band prop¬ 
erties p!2lll3[ll4lll5lll6j . Other studies for unsaturated 
granular media observe fluid depletion in shear bands 
mM- However, there is no theoretical framework or 
concrete model available yet that defines the exact cor¬ 
relation between the micro parameters like the liquid 
bridge volume and the surface tension of the liquid with 
the steady state cohesion. 

In order to develop a micro-macro correlation for 
the liquid bridge contact model, we initially study the 
structure of the micro contact model. How is the struc¬ 
ture of the liquid bridge contact model affected by the 
microscopic parameters? How does this influence the 
steady state cohesion? Here we study in detail on the 
effect of these parameters on the macro results. For ex¬ 
ample, the effect of maximum interaction distance, or 
the distance at which the liquid bridge between two in¬ 
teracting particles ruptures, is studied by varying the 
liquid content. On the other hand other parameters like 
surface tension of the liquid and contact angle affect the 
magnitude of force acting between the particles when in 
contact naiE]. Various surface tension of liquids give 
a large scale variation of the capillary force and this 
allows us to study the effect of maximum force on the 
macroscopic properties. Furthermore, in the consecu¬ 
tive analysis, we re-obtain the macro-rheology results 
in the shear band center from the torque, torque being 
an experimentally measurable quantity. 

The liquid bridge interactions between the particles 
are defined by the free-surface equilibrium shapes and 
stability of the bridge configuration between them [20l 
mm- Phenomenologically, even the simplified models 
of liquid bridges are quite complex in nature. In order 
to improve the computational efficiency for wet granu¬ 
lar materials, we replace the non-linear interactions of 
liquid bridges with a simpler linear one. But in what 
way can a non-linear model like the liquid bridge con¬ 
tact model be replaced by a linear model? When can 
we say that the two different contact models are analo¬ 
gous? Therefore, we compare the realistic liquid bridge 
model with an equivalent simple linear irreversible con¬ 
tact model [23] that would give the same macroscopic 
effect. 


The results in this paper are organized in three main 
parts. In Sec. |3.1| of this paper we study the effect 
of varying liquid bridge volume and surface tension of 
the liquid on the macroscopic properties, the focus be¬ 
ing to find a micro-macro correlation from this study. 
Most strikingly, we see a well defined relationship be¬ 
tween these micro parameters and the macro properties 
like the steady state cohesion of the bulk material and 
macro-torque required under shear, neglecting the effect 
of fluid depletion in shear bands [num in quasistatic 
flow. In Sec. |3.2| of this paper we show the derivation of 
macro torque from the boundary shear stress. In this 
section we also compare this torque with the torque cal¬ 
culated from forces due to contacts on the wall particles. 
In Sec. iof this paper, we discuss about the analogy of 
two different contact models, with a goal to understand 
which parameters at microscopic scale would give the 
same macroscopic behavior of the system. 


2 Model System 

2.1 Geometry 

Split- Bottom Ring Shear Cell: The set-up used for sim¬ 
ulations consists of a shear cell with annular geometry 
and a split in the bottom plate, as shown in figure 
Some of the earlier studies in similar rotating set-up in¬ 
clude [2411251126] . The geometry of the system consists 
of an outer cylinder (radius Rq = 110 mm) rotating 
around a fixed inner cylinder (radius Ri = 14.7 mm) 
with a rotation frequency of frot = 0.01 s“^. The gran¬ 
ular material is confined by gravity between the two 
concentric cylinders, the bottom plate, and a free top 
surface. The bottom plate is split at radius Rs = 85 
mm into a moving outer part and a static inner part. 
Due to the split at the bottom, a shear band is formed 
at the bottom. It moves inwards and widens as it goes 
up, due to the geometry. This set-up thus features a 
wide shear band away from the wall, free from bound¬ 
ary effects, since an intermediate filling height {H = 40 
mm) is chosen, so that the shear band does not reach 
the inner wall at the free surface. 

In earlier studies similar simulations were 

done using a quarter of the system (0° < 0 < 90°) 
using periodic boundary conditions. In order to save 
computation time, here we simulate only a smaller sec¬ 
tion of the system (0° < 0 < 30°) with appropriate pe¬ 
riodic boundary conditions in the angular coordinate, 
unless specified otherwise. We have observed no notice¬ 
able effect on the macroscopic behavior in comparisons 
between simulations done with a smaller (30°) and a 
larger (90°) opening angle. Note that for very strong 
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Fig. 1 Shear cell set-up. 

attractive forces, the above statement is not true any¬ 
more. 


2.2 Microscopic model parameters 

In presence of a small amount of liquid in a dense gran¬ 
ular material, bridges are formed at the contact points 
between the particles. The surface energy of these bridges 
leads to an attractive force between the particles, which 
is absent in dry granular materials. Thus, wetting changes 
a granular system from one with only repulsive inter¬ 
particle interactions to one with both repulsive and at¬ 
tractive interactions [29]. With the change in micro¬ 
scopic physical interactions in wet granular materials, 
the macroscopic behavior is also expected to differ from 
the dry materials. Therefore, we choose to vary some of 
the characteristic specifications of a liquid bridge model 
to understand the effect on macroscopic properties. All 
the particle specifications and the fixed interaction pa¬ 
rameters for the contact models are given in table 
All the variable interaction parameters which include 
the liquid bridge volume Vb and the surface tension of 
the liquid 7 are discussed in this section. 

2.2.1 Bulk saturation and liquid bridge volume 

The bulk material can be characterized by different 
states such as the dry bulk, adsorption layers, pendu¬ 
lar state, funicular state, capillary state or suspension 
depending on the level of saturation [301131] . In this pa¬ 
per we intend to study the phenomenology of liquid 


bridge between particles in the pendular state, where 
the well separated liquid bridges exist between particle 
pairs without geometrical overlap. In this section, we 
discuss about the critical bulk saturation of granular 
materials and the corresponding liquid bridge volumes 
in the pendular state. 

The bulk saturation S'* is defined as the ratio of 
liquid volume to void volume of the bulk [3^1331134] . 
The demarcation between the pendular state and the 
more saturated funicular state is given by the saturation 
S* « 0.3 [32]. For each particle pair with a liquid bridge, 
a dimensionless volume (^* can be defined as the ratio 
of the volume of the liquid bridge at the contact, V 5 lo 
the volume of the two contacting particles, 2V^: 


^ 21 /p 2(fd3) 

Assuming the liquid is homogeneously distributed through¬ 
out the material, the bulk saturation S* is obtained 
from the dimensionless volume (/:?* and the bulk poros¬ 
ity e from the following equation [3^1331134] : 


5 * = 2-^^* (2) 

With a bulk porosity of the material e = 0.4 and a mean 
particle diameter dp of 2.20 mm, the maximum liquid 
bridge volume in the pendular regime is approximately 
284 nl. In order to study the influence of liquid content 
on the macroscopic properties, we analyzed the system 
for the following set of liquid bridge volumes W- 


Vb e {0, 1, 2, 4.2, 8, 14, 20, 75, 140, 200} nl, (3) 


which are seen to be well within the pendular regime. 
We also calculate the liquid volume as a percentage of 
the total volume of the system (Vt) based on the number 
of contacts. The number of contacts represented as Cl 
increases with increasing liquid bridge volume in the 
system and is measured approximately: 


Cl G {33010 , 36214, 36855, 37585, 38306, 39101, 
39511, 41526, 42595, 43328}, (4) 

Therefore, the volume percentage of liquid in the sys¬ 
tem is given by cpb = and is approximately equal 

to: 

ipb e {0 , 0.03, 0.07, 0.15, 0.29, 0.52, 

0.75, 2.94, 5.63, 8.18}. (5) 

In order to investigate the functional form of steady 
state cohesion beyond this state, a few more simulations 
for higher Vb are done: 

Vb e {500, 1000} nl, (6) 

for which the pendular assumption is not valid any¬ 
more. 
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Table 1 Model parameters 


Parameter Symbol Value 


Sliding friction coefficient 

Elastic stiffness 

Viscous damping coefficient 

Angular frequency 

Particle density 

Mean particle diameter 

Contact angle 



0.01 

k 

120 Nm-i 

7o 

0.5x10“^ kgs 

UJ 

0.01 s-i 

P 

2000 kgm-3 

dp 

2.2 mm 

0 

20° 


2.2.2 Surface tension of liquid 

Surface tension results from the greater attraction of 
liquid molecules towards each other than towards air. 
It is the tendency of liquids to lower their state of en¬ 
ergy which makes it acquire the least possible surface 
area at the surface with higher inter liquid molecules 
attraction. As a result, cohesive properties of liquids 
are reflected in surface tension which makes it an inter¬ 
esting parameter to study. This effect will be discussed 
in detail in Sec. inn The effect of surface tension on 
the macroscopic properties is studied for the following 
range of surface tension values: 

7 e {0, 0.020, 0.040, 0.060} Nm-r (7) 

Surface tension of most of available liquid-air interfaces 
at 20°C are in this range. To investigate the functional 
behavior of steady state cohesion beyond this state, a 
few more simulations for higher 7 are done: 

7 e {0.01, 0 . 10 , 0.50, 1.00} Nm“t ( 8 ) 


2.3 Liquid bridge contact model 


The contact and non-contact forces for interacting par¬ 
ticles can be described by a combination of an elastic 
contact model for the normal repulsive force and a non¬ 
linear irreversible adhesive non-contact model for the 
adhesive force. Figure [^represents a sketch of the com¬ 
bined liquid bridge contact model as a function of the 
overlap between the two particles. The liquid bridge ad¬ 
hesive force acts between the particles once the contact 
is established. When the particles are in contact, the at¬ 


tractive force is given by Eq. (13). This is independent 


of the liquid bridge volume and depends on the surface 
tension of the liquid, radius of particles and contact 
angle. There is no cohesive force between the particles 
during approach. As the liquid bridge only forms once 
the particles come in contact with each other, the co¬ 
hesive force starts acting and remains constant during 
overlap between particles > 0. Normal contact re¬ 
pulsive force acts between the particles in contact in 



Fig. 2 Liquid capillary bridge model. The red lines represent 
the loading direction, the blue line represents the unloading 
direction when the particles are in contact and the brown line 
represents the unloading for the non-contact particles with 
short-range interaction force. 


addition, given by: 

fn = kS + 7o(5, (9) 


where k is the elastic stiffness, 7 o is the viscous damping 
coefficient and 6 is the overlap between the particles. 
The normal contact forces for the liquid bridge model 


are explained in Sec. 2.3.1 


2.3.1 Liquid bridge capillary force model 

The capillary pressure difference sustained across the 
liquid-air interface due to surface tension can be de¬ 
scribed by the non-linear Laplace-Young equation [22]. 
This relates the pressure difference to the shape of the 
surface under the criterion of minimum Gibbs free en¬ 
ergy [35] . The capillary force in a pendular bridge orig¬ 
inates from the axial component of this force. Another 
component that contributes to the capillary force is due 
to the hydrostatic pressure. Many previous studies have 
calculated capillary forces based on the numerical solu¬ 
tion of the Laplace-Young equation and also reported 
experimental results nana. The magnitude of liquid 
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bridge capillary force depends on the volume of the liq¬ 
uid bridge between the particles, the contact angle 6 >, 
surface tension 7 , the effective radius of the particles 
r and the separation distance S, S = —S. With these 
parameters we approximate the inter-particle force fc 
of the capillary bridge according to m- The experi¬ 
mental results are fitted by a polynomial to obtain the 
dependence of capillary forces on the scaled separation 
distance. During approach of the particles, the normal 
contact force for this model is given by: 


/ = 


( 


0 



if <5 < 0; 

fn if (^ > 0. 


( 10 ) 


During separation of the particles, the normal con¬ 
tact force for this model is given by; 


/ 


'0 if(5<-S'c; 

< -fa if -Sc < (5 < 0; 


( 11 ) 


where /„ is the normal repulsive force given by Eq. 
The adhesive force for the liquid bridge model is 
the capillary force given by: 


fa 


if a) 


liq 


( n max\ ( 2r\ 

Ua )liq\dp) 

1 +LOSS'+ 2.55'2’ 


( 12 ) 


where the separation distance is normalised as S = 
S^/{r]V^^ S being the separation distance. The max¬ 
imum capillary force between the particles when they 
are in contact (S = 0 ) is given by: 




(13) 


where, dp is the mean particle diameter. The effective 
radius of two interacting spherical particles of different 
size can be estimated as the harmonic mean of the two 
particle radii according to the Derjaguin approximation 
[36], yielding the effective radius: 


2rirj 


(14) 


however, the mean size is not varied here. This model 
equation is applicable for mono-disperse particles [ 12 ] 
HD which has been actually extended to poly-disperse 
system of particles Ref. mi. As proposed by the 
critical separation distance Sc between the particles be¬ 
fore the bridge ruptures is given by: 


(5c)a,= (l + 0L^" 


(15) 


The liquid bridge capillary force as a function of sep¬ 
aration distance is shown in figure for different liq¬ 
uid bridge volumes. The liquid bridge capillary force 


- 0.2 
-0.4 
- 0.6 
- 0.8 
-1 
- 1.2 
-1.4 
- 1.6 

-1.8'-'-'-'-^^^-' 

-0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 

6 * 



Fig. 3 /e* as a function of (5*. Different colors represent 
different liquid bridge volumes. 


decreases in magnitude with increasing separation dis¬ 
tance between the particles till the bridge ruptures. The 
rupture distance is proportional to ' as stated in Eq. 

(itl). 


2.3.2 Linear irreversible contact model 



Fig. 4 Linear irreversible contact model. The red lines rep¬ 
resent the loading direction, the blue line represents the un¬ 
loading direction when the particles are in contact and the 
brown line represents the unloading for the non-contact par¬ 
ticles with short-range interaction force. 


In Sec. we introduce a simple linear irreversible 
contact model as proposed by [23] and shown in fig¬ 
ure which can be compared with the non-linear liquid 
bridge interaction model. For the linear irreversible con¬ 
tact model, the normal forces between particles during 
approach and separation are given by Eqs. (10) and 
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(11) respectively, where for the linear irreversible con¬ 
tact model, 


3.1 Steady state cohesion and its correlation with 
liquid bridge volume and surface tension 


fa 


(/aW = (/a'"“")a„ + M, 


(16) 


{Schn = ifa^^niinlkc, (17) 

Where i/mun is the maximum adhesive force and 
kc is the adhesive stiffness. The tangential force contact 
model is explained in details in our earlier studies EH- 


2.4 Dimensional analysis 

To formulate all the modeling equations in a construc¬ 
tive way, we express them in nondimensionalized form. 
All the length scale parameters are scaled by the mean 
particle diameter dp = 2.20 mm. The forces are scaled 
in terms of the gravitational force acting on a single 
particle fg = Vppg ^ 1.0939 x 10“^ N. Table I shows 
all the parameters in their dimensionless form and the 
corresponding scaling terms used in the equations. The 
angular rotation of the shear cell after a given time 
to study the dynamic evolution of torque is scaled in 
terms of radians covered in one complete rotation (27r). 
The dynamics of the system can be characterized by 
the time scale defined by the contact duration between 
two particles tc = where rrip is the mean mass 

of a particle. Since we do all our macro-rheology analy¬ 
sis in steady state, characterization of dynamics of the 
system is not required. The main objectives of nondi- 
mensionalization is to simplify the equations in terms of 
unit less quantities and define the system intrinsically. 


3 Micro macro transition 

To extract the macroscopic properties, we use the spa¬ 
tial coarse-graining approach detailed in [3811391140] . The 
averaging is performed over toroidal volume, over many 
snapshots of time assuming rotational invariance in the 
tangential ^-direction. The averaging procedure for a 
three dimensional system is explained in [3^140] . This 
spatial coarse-graining method was used earlier in [Tll^ 
EillSHlllo]. The simulation is run for 200 s and temporal 
averaging is done when the flow is in steady state, be¬ 
tween 80 s to 200 s, thereby disregarding the transient 
behavior at the onset of the shear. 


In earlier studies the shear band region 

was identified by the criterion of large strain rate, e.g. 
higher than a critical strain rate of 0.08 s“^. In this 
paper, the shear band center region is defined by strain 
rates higher 80% of the maximum for different heights 
in the shear cell. Figure displays the dependence of 
scaled yield stress r* for the particles in the shear band 
region on scaled pressure P* for 75 nl liquid bridge vol¬ 
ume. A linear trend is observed neglecting the different 
behavior for data at very low pressure (P* < 4.42). 
This is fitted well by a linear function: 

r*=/iP*-hc* (18) 


where fi is the macroscopic friction coefficient and c* is 
the steady state cohesion obtained from the plot. Next, 
we fit the data for shear stress as a function of pressure 
as given by Eq. (18) and obtain the value of steady state 
cohesion and macroscopic friction g. The macroscopic 
friction coefficient is constant for lower surface tension, 
including 7* = 0 for linear elastic model (not shown 
in figure), but increases for 7* > 2 for a given liquid 
bridge volume as shown in figure When the surface 
tension of the material is very high (7* > 1.00), ma¬ 
terials protrude out of the top surface to form a hump 
in the region of the shear band (data not shown). For 
our analysis of surface tension in the range 0.020-0.040 
Nm“^, the macroscopic friction coefficient is constant 
8it g 0.15. In this range, the macroscopic friction co¬ 
efficient is also independent of the liquid bridge volume 
as shown in figure 

For dry cohesionless systems, the dependence of shear 
stress on pressure is linear without an offset, i.e. c* = 
0. In the presence of interstitial liquid between the par¬ 
ticles in the pendular regime, cohesive forces increase 
with increasing liquid bridge volume. This results in a 
positive steady state cohesion c* as given by Eq. (18), 
see figure 

Earlier studies on wet granular materials have shown 
that the presence of liquid bridges between the parti¬ 
cles results in an increasing steady state cohesion of 
the materials [T^[T3l[T^I^ . Our earlier studies show 
that the steady state cohesion c* increases non-linearly 
with increasing liquid bridge volume. Here, the steady 
state cohesion is studied in more detail, including very 
small liquid bridge volumes, including the (practically 
impossible) limit of 0 nl liquid bridge volume as given in 
Eq. (§. Note that there is a finite cohesive strength for 
W ^ 0 nl liquid bridge volume. This is due to the micro¬ 
scopic capillary bridge force that acts between particles 
even at 0 nl liquid bridge volume as given by Eq. (13). 
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Table 2 Non-dimensionalization of parameters 


Parameter 

Symbol 

Scaled term 

Scaling term 

Capillary force 

fc 

fc^ 

f9 

Particle overlap 

S 

(5* 

dp 

Shear stress 

T 

r* 

fgldp^ 

Pressure 

P 

P* 

fgjdp^ 

Steady state cohesion 

C 

c* 

fg/dp^ 

Liquid bridge volume 

Vb 

H* 

dp^ 

Surface tension 

7 

7 * 

fg/dp 

Rupture distance 

5e 

5c* 

dp 

Torque 

T, 

rp * 

fgdp 

Angular rotation 

^vot 

0 * 

^rot 

27r 

Adhesive Energy 

E 


fgdp 




Fig. 5 Shear stress r* plotted against pressure P*. The dot¬ 
ted line represents the fitting function as given by Eq. ( |18| ) 
for P* > 4.42 Pa where /r = 0.15 is the macroscopic friction 
coefficient, c* = 0.2655 for Vb = 75 nl and 7 = 0.020 Nm“^. 


0.35 

0.3 


0.25 

=L 

0.2 


0.15 


10’ 10“ 10’ 10^ 

T 

Fig. 6 Macroscopic friction coefficient /r as a function of 7 * 
for Vb = 75 nl. The solid symbols represent the range of sur¬ 
face tension for our simulations below. 



Fig. 7 Macroscopic friction coefficient // as a function of 5c* 
for 7 = 0.020 Nm“^. 


This is called the steady state critical cohesion cq* for 
a given surface tension of liquid. This value depends on 
the maximum force acting between two particles when 


they are in contact as given by Eq. (13). The additional 
cohesion for higher liquid bridge volume is due to the 
non-contact capillary forces between the particles that 
are active upto the distance when the liquid bridge rup¬ 
tures. This is dependent on the surface tension of the 
liquid and the volume of the liquid bridge. Thus, the 
steady state cohesion of granular materials for a given 
liquid bridge volume can be written as: 


Cq + C 


y* 


(19) 


where c'* is the additional cohesion for liquid bridge 
volume Vb > 0. Figure a) shows (c* — cq*) as a linear 
function of 5c*, fitted by: 


c* - Co* = a5c* (20) 

where a = 0.9805 for 7 = 0.020 Nm“^. In the next 
section we study the dependence of this constant on 
the surface tension of liquid. 
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S " 


c 



Fig. 8 (a) c* — Co* as a function of Sc* for 7 = 0.020 Nm“^. 
The dotted line represents the fitting function given by Eq. 
( |20| ). The data with solid symbols represent the liquid bridge 
volume outside the pendular regime, (b) c* as a function of 7 * 
for Vb = 75 nl . Th e dotted line represents the fitting function 
given by Eq. ( [2T| . 


in terms of maximum interaction distance ^Sc* between 
the particles and the maximum adhesive force expressed 
in terms of surface tension of the liquid 7 *. So in the 
later sections of this paper we study the dependence of 
macroscopic parameters on the micro parameters S'c* 
representing scaled rupture distance and 7 * represent¬ 
ing scaled maximum force for all contact models. 

Figure shows the dependence of ^ on S'c* for 
different surface tension of liquid. The scaled steady 
state cohesion is a linearly dependent on the rupture 
distance as shown in the figure. This can be fitted by a 
straight line equation given by: 


c* - Co _ c* - Co _ 

7* “ ifa^nu^/i^COsO) - 


(23) 


where p = 2.1977 as obtained from the fitting shown in 
figure the offset is very small and can be neglected. 
This subsection shows that the macroscopic character- 



Figure [^b) shows the dependence of steady state 
cohesion on 7 * for Vb = 75 nl. The steady state cohesion 
can be described by: 

Inc* = (tin 7 * + A: ( 21 ) 

where a ~ 1.00, k = —0.4240. Therefore, the steady 
state cohesion is linearly proportional to the surface 
tension and can be written as: 


= 67 * 


( 22 ) 


where b = exp{k). The above equation is valid in the 
limit of zero surface tension ( 7 * = 0 ) which represents 
the simple linear elastic contact model. For higher sur¬ 
face tension of liquid, the results deviate from the fitted 
function of linear dependence as seen from figure [^b). 
As given by Eq. (20) and ( 21 ), the steady state co¬ 
hesion is dependent on liquid bridge volume expressed 


Fig. 9 as a function of Sc* for different surface ten¬ 

sion of liquid. The dotted line represents the fitting function 
given by Eq. ( [^ . 


istics of the liquid bridge model are determined by the 
maximum interacting force between the particles and 
the rupture distance. The steady state cohesion scales 
linearly with the surface tension of liquid i.e. the maxi¬ 
mum force between the particles. For a given maximum 
force, the cohesion scaled with the surface tension of liq¬ 
uid is also a linear function of the rupture distance of 
the liquid bridge. 


3.2 Macroscopic torque analysis from the microscopic 
parameters 

The strength, cohesion and flow properties of granu¬ 
lar materials are strongly influenced by the presence 
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of capillary cohesion. Due to the cohesive properties 
of these wet materials, the shear stress increases and, 
as a result, partially saturated wet materials require 
higher torques for deformation (shear) e.g. in a shear 
cell. Loosely speaking, torque is a measure of the shear 
stress or force acting on the particles at the wall and 
thus can be used to find an estimate of shear stress in 
the shear band. To study solely the effect of capillary 
cohesion on the torque, the other parameters like the 
particle friction is kept very small in our simulations, 
with jjLp = 0.01. Earlier studies p^l2 71141 [[42] show that 
the average torque acting on the rotating part of the 
shear cell increases with increasing moisture content. 
In this section we perform a detailed analysis of the 
macroscopic torque as a function of the micro param¬ 
eters in order to understand its connection with the 
steady state cohesion of the material. 

The walls and the bottom plates of the shear cell 
consist of particles with a prescribed position. The par¬ 
ticles forming the inner wall are stationary while the 
particles forming the outer wall rotate around the z-axis 
with frequency frot- All the particles forming the inner 
and outer wall are identified as Cinner and Couter, respec¬ 
tively. The macroscopic torque is calculated based on 
the contact forces on the fixed particles on the moving 
(outer) and stationary (inner) parts of the shear cell. 
Thus the net inner and outer torque are calculated by 
summing up the torques for all the contacts with re¬ 
spect to the axis of rotation of the shear cell. The net 
torque is obtained from the difference between the outer 
wall torque and the inner wall torque. We multiply the 
total torque by a factor of in order to get the torque 
for the whole system from the obtained torque of our 
simulations in a 30° section. Thus the torque is given 
by: 


27r 


N 


T= Cij xfij)- 

LVZ-.Z^^eCouter 
' 1=1 


N 


(EE. 


2=1 


'ieCi. 


xfij)], (24) 


where N represents the number of particles, Cij is the 
position of the contact point and ^ij is the interaction 
force. Only the ^(-component of the torque vector (T^) 
is of interest as required for shearing the cell in angular 
direction. 

We compare our results with the experimental re¬ 
sults as given by [43] from the evolution of torque as 
a function of the angular rotation as shown in figure 
[Tq| This is in good agreement with the magnitude and 
angular rotation required for steady state torque evolu¬ 
tion as given in m, considering the different rotation 
rate and different friction in the systems. 



Fig. 10 T^* as a function of scaled angular rotation Orot* 
for surface tension of liquid 7 = 0.020 Nm“^ for V 5 = 4.2 nl 
(blue) and V 5 = 200 nl (red). 


Figure [TT] shows T^* as a function of 7* for different 
liquid bridge volumes. We observe that the resultant 
torque depends linearly in the surface tension of the 
liquid. The fit parameter I from the figure, the rate of 
increase of torque with surface tension, depends on the 
liquid bridge volume. 

Next, we compare the results of the steady state 
cohesion as obtained from the fitting function explained 
in Sec. 3.1 with the calculated (measured) torque. We 
write the scalar form of the torque on the wall derived 
from steady state cohesion as 


rji macro _ 

/ r dA — r dA 


JAo JAi 


{f^Pavg c): 


(25) 


where Ao denotes the outer wall surface, Ai denotes the 
inner wall surface and Pavg is the mean pressure inside 
the shear band approximately 250 Pa for a filling height 
of 39 mm. Eq. ( [2^ can be simplified to the form: 


= + (26) 

where M = [27rH{Ro^—Ri‘^)-\-^7r{Ro^-\-Ri^— 2Rs^)] ~ 
0.0031 m^ for the given geometry is equal to fitting pa¬ 
rameter t/{jiiPavg)^ t is the fit parameter, see figure [TT| 
Assuming Tz = ^ an equivalent steady state co¬ 

hesion as obtained from the calculated torque can be 
given as: 


Ceq = Tz/M - jlPavg {27) 

Figure [T^ shows the dependence of the non-dimen- 
sionalised value on 5'c* for different surface 

tension, (ceg)o* is the equivalent steady state cohesion 
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as obtained from Eq. (27) for the torque of a 0 nl liquid 


bridge. This can be fitted by a straight line: 


Ceq* - (Ce,)o* 

'-y* 


= e5c*, 


(28) 


where e = 2.0062 is a fit parameter, see figure [T^ and 


the offset is very small and can be neglected. Eq. (28) 


shows equivalent steady state cohesion as obtained from 
the torque is also linearly dependent on Sc ^ The fitting 
parameter e of this equation shows a close similarity 


with the fitting parameter p of Eq. (23). Alternatively. 


figure [T^ shows a comparison of the two torque given by 
the scalar ^-component of Eq. (24) and Eq. (26) for sur¬ 


face tension of liquid 0.020 Nm These results show 
that the steady state cohesion and torque are related 


by Eq. (26). 


In conclusion, this subsection shows that the mea¬ 
sured torque can be translated to the local steady state 
macro-rheology parameters via a simple factor M (a 
measure of the resultant arm-length times surface area) 
which depends only on the geometry of the system. 




Fig. 12 ^as a function of for different surface 
tension of liquid where Ceq is given by Eq. ( |27| ). The dotted 
line represents the fitting function as given by Eq. ( |28| ). 



Fig. 13 Torque calculated numerically scaled as T;^* as com¬ 
pared with the scalar form of scaled macro torque 
as calculated from the wall shear stress as given by Eq. ( |26| ). 


Fig. 11 T,* as a function of 7 *. The dotted lines represent 
the fitting functions for different liquid bridge volumes given 
by equation Tz* = Ij* + t where t = 4.964 xlO^ and I in¬ 
creases with increasing liquid bridge volume. 


4 An analogous linear irreversible contact 
model for cohesive particles 

In this section we aim to determine the key microscopic 
parameters for a linear irreversible contact model [23] 
that is macroscopically analogous to the liquid bridge 
contact model used before. An explanation of the lin¬ 
ear irreversible contact model is given in [23|. Unlike the 


liquid bridge contact model, the force for the linear irre¬ 
versible contact model is simple and faster to compute. 
Figure shows the force-overlap distribution for the 
two contact models showing the loading and unloading 
directions of forces which are reversible at ^* > 0 and 
irreversible at (5* < 0. 

As discussed in Sec. |3.1[ the steady state cohesion 
for the liquid bridge model is controlled by the rup¬ 
ture distance of the liquid bridge, which is proportional 
to the liquid bridge volume, and the magnitude of the 
maximum interaction force, which is governed by the 
surface tension of the liquid. Assuming that the non¬ 
linear liquid bridge capillary force can be replaced by 
a simple irreversible linear adhesive force between the 
particles with the same macro characteristics, we com- 
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Fig. 14 Force-overlap diagram for the liquid bridge model 
(blue) as compared with the linear-irreversible contact model 
(red). The arrow shows the loading and the unloading direc¬ 
tions for all forces. The schematic diagram for the same are 
given in figures and respectively. 

pare the steady state cohesion of the two models in Sec. 

HU 


4.1 Equal maximum force and interaction distance 

The key parameters that define the cohesive force of 
a linear irreversible contact model are the maximum 
adhesive force and the adhesive stiffness, see Eq. (16). 
Several simulations have been run for the linear irre¬ 
versible contact model in the same numerical set-up 
with the same maximum adhesive force as used in the 
liquid bridge model and ad- 

hesive stiffness that would result in the same interaction 
range for different liquid bridge volumes for different 
surface tension of liquid. The force-overlap for contacts 
with d* < 0 for the two comparable contact models 
with equal interaction distance are shown in figure 
The adhesive stiffnesses that are equivalent to the liquid 
bridge volumes as given by Eq. for surface tension 
7 = 0.020 Nm“^ for equal interaction distance are given 
by; 

ka e {0.21 , 0.26, 0.41, 0.46, 

0.56, 0.69, 0.88, 1.11, 00} Nm“^ (29) 

The results for the steady state cohesion c*, as scaled 
by 7* for the liquid bridge model and the linear irre¬ 
versible model are shown in figure The results are 
not really analogous as seen from the figure as the inter¬ 
cepts for the fitting lines of the two models are different, 
while they are parallel. The fitting parameters for the 



5* 


Fig. 15 Scaled adhesive force /* (/a* for linear adhesive 
model and /c* for liquid bridge model) as a function of (5* 
for the linear irreversible contact model (red), compared with 
the liquid bridge model (blue), for equal maximum force and 
equal interaction distance. The yellow line represents the force 
for the liquid bridge contact model for mean particle diameter 
dp as a function of (5*. The arrow shows the loading and the 
unloading directions for the short-range forces. 

relation: 
c* — Cn* 

—^ = gSc* + h (30) 

are g = 2.1716 and h 0 for the liquid bridge con¬ 
tact model, g = 2.0984 and h = 0.2226 for the linear 
irreversible contact model. 

So for a given liquid bridge volume and a given 
surface tension of liquid, the linear irreversible contact 
model with the same maximum force and same inter¬ 
action distance has a higher cohesion. 

4.2 Equal maximum force and adhesive energy 

Equal maximum force and interaction distance was dis¬ 
cussed in Sec. 4.1, but here the steady state cohesion 
for the two models with an equal maximum adhesive 
force and equal adhesive energy E* are considered. The 
adhesive energy for a given contact model is obtained 
by the total area under the force-overlap distribution, 
see figure A linear model analogous to the liquid 
bridge contact model is obtained with the equal max¬ 
imum force with surface tension 7 = 0.020 Nm“^ and 
the adhesive stiffness adjusted to have the equal adhe¬ 
sive energy: 

ka G {0.25, 0.29 , 0.39, 0.74, 0.84, 

1.10, 1.49, 2.11, 2.95, 00} Nm“^ (31) 















12 


Sudeshna Roy et al. 



6 * 

Fig. 16 Scaled adhesive force /* (/a* for linear adhesive 
model and /c* for liquid bridge model) as a function of (5* 
for linear the irreversible contact model (red), compared with 
the liquid bridge model (blue), for equal maximum force and 
equal adhesive energy dissipated per contact. The yellow line 
represents the force for the liquid bridge contact model for 
mean particle diameter dp as a function of (5*. The arrow 
shows the loading and the unloading directions for the short- 
range forces. 


equal energy has a lower interaction distance. The func¬ 
tional behavior of the steady state cohesion using the 
linear irreversible contact model for small interaction 
range can be understood from this. As observed from 
figure the cohesion is a non-linearly dependent on 
the rupture distance 5'c* at low interaction distance and 
becomes linear for higher range. 

Figure shows the dependence of steady state co¬ 
hesion on total adhesive energy for the liquid bridge 
model, compared with the two cases of linear irreversible 
contact model with equal interaction distance (red) and 
equal adhesive energy dissipated per contact (green). 
As seen from the figure, for a given maximum force 
which is determined by the surface tension of the liq¬ 
uid, the steady state cohesion c* is equal for the the 
liquid bridge model and the linear irreversible model 
with equal energy. The steady state cohesion for the lin¬ 
ear irreversible model with equal interaction distance is 
higher as it has higher adhesive energy than the liquid 
bridge model. However, all the data for the three cases 
as explained above collapse and functionally behave the 
same. 



Fig. 17 as function of Sc * for the liquid bridge model 

(blue) and the linear irreversible model with equal interaction 
distance (red) and equal adhesive energy dissipated per con¬ 
tact (green) for 7 = 0.020 Nm“^. The dotted a nd the solid 
lines represent the fitting function given by Eq. ( |30| ). 


Fig. 18 c* as a function of E* for the liquid bridge model 
(blue) and the linear irreversible model with equal interac¬ 
tion distance (red) and equal adhesive energy dissipated per 
contact (green) for 7 = 0.020 Nm“^. 


The force-overlap for contacts with (5* < 0 for the 
two comparable contact models with equal adhesive en¬ 
ergy are shown in figure pT) Figure shows the depen¬ 
dence of on rupture distance Sc* for the liquid 

bridge model (blue), compared with the two cases of 
the linear irreversible contact model with equal interac¬ 
tion distance (red) and equal adhesive energy dissipated 
per contact (green). The linear irreversible model with 


4.3 Different maximum force for the two contact 
models 

In the earlier subsections, results show that for a given 
maximum force the steady state cohesion for the two 
contact models functionally behave the same under equal 
force and equal energy conditions. To study the func¬ 
tional form for the two models under different maxi- 
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mum force conditions, we compare the macroscopic be¬ 
havior of the linear model to the liquid bridge model 
results for different surface tension. Linear model sim¬ 
ulations equivalent to surface tension 0.040 Nm“^ and 
0.060 Nm“^ are run with an equivalent adhesive stiff¬ 
ness 2 times and 3 times of that given by Eq. keep¬ 
ing the interaction distance the same. Figure^] shows 
a comparison of the force-overlap for the two contact 
models for surface tension of liquid 0.020 Nm“^ and 
0.040 Nm-h 





8 * 


Fig. 19 Scaled adhesive force /* (/a* for linear adhesive 
model and fc"" for liquid bridge model) as a function of 8* 
for the linear irreversible contact model (red), compared with 
the liquid bridge model (blue) for different maximum force 
and equal interaction distance. The yellow lines represent the 
force for the liquid bridge contact model for mean particle 
diameter dp as a function of 8*. The arrow shows the loading 
and the unloading directions for the short-range forces. 

Figure [2Q|^a) shows the dependence of steady state 
cohesion on the adhesive energy dissipated by the parti¬ 
cles per contact for different for the liquid bridge 

model and the linear model. For the same energy dissi¬ 
pated per contact, a higher surface tension of the liquid 
results in a higher macroscopic cohesion. Figure [2Q|^b) 
shows that ^ is a function of ^ for a given surface 
tension, or maximum force. 

5 Conclusion 

We observed a correlation between the steady state 
cohesion and the microscopic parameters of the liq¬ 
uid bridge model. The micro-parameters are the liquid 
bridge volume, the liquid surface tension, the contact 
angle (which was kept constant) and the size of particles 
(i.e. curvature, which was also not varied). A detailed 
study of the effect of liquid bridge volume and surface 




Fig. 20 (a) c* as a function of F* for different surface ten¬ 
sion of liquid, (b)^ as a function of ^ for different sur¬ 
face tension of liquid as compared with the linear irreversible 
model. Different symbols denote o liquid bridge model and V 
linear irreversible model. 


tension of liquid is done in this paper. These micro¬ 
scopic parameters control the cohesion in wet granular 
materials in different ways. We found that the steady 
state cohesion of the system is proportional to the maxi¬ 
mum adhesive force, which was varied by modifying the 
surface tension. On the other hand, we found that the 
steady state cohesion is also linearly dependent on the 
maximum interacting distance between the particles, 
which depends on the volume of the liquid bridge. From 
these results we have obtained a good micro-macro cor¬ 
relation between the steady state cohesion and the mi¬ 
croscopic parameters studied. 

We analyzed the effect of cohesion on the wall torque 
required to rotate the system at a given rate. The torque 
and the steady state cohesion of the system are pro¬ 
portional and show similar linear dependence on the 
microscopic parameters. 
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Finally, an analogy was established between the liq¬ 
uid bridge model and a simpler linear irreversible con¬ 
tact model; even though these two models have different 
micro-macro correlations, the steady state cohesion for 
the two models are the same if the maximum force and 
the total adhesive energy dissipated per contact for the 
two models are matched, irrespective of the shape of 
the attractive force function acting between the parti¬ 
cles. In this way one can always replace a non-linear 
liquid bridge force by a simpler, faster to compute lin¬ 
ear one and obtain identical macroscopic properties in 
less computational time. 

Furthermore, results for the two types of contact 
models with equal energy and different magnitude of 
maximum force show that they have different steady 
state cohesion. Therefore, adhesive energy is not the 
sole microscopic condition for the two contact models 
to have same steady state cohesion. Instead, we observe 
that both adhesive energy and cohesion scale linearly 
with the maximum adhesive force. In this way we can 
determine the steady state cohesion from the two mi¬ 
croscopic parameters, the adhesive energy and the max¬ 
imum force. 

In this paper our study was focused on the micro¬ 
macro correlations and comparing different contact mod¬ 
els. It would be interesting to also study the forces and 
their probability distributions for wet cohesive systems. 
Future studies will aim at understanding the micro¬ 
scopic origin and dynamics of the contacts and liquid 
bridges throughout the force network(s) and also the 
directional statistics of the inter-particle forces inside a 
shear band. The effect of liquid migration on the macro 
properties and a continuum description for a similar 
model will be studied in the near future. 

Acknowledgements We acknowledge our financial support 
through STW project 12272 ’’Hydrodynamic theory of wet 
particle systems: Modeling, simulation and validation based 
on microscopic and macroscopic description.” 


References 

1 . A. Singh, V. Magnanimo, K. Saitoh, S. Lading, Physical 
Review E 90(2), 022202 (2014) 

2. N. Mitarai, H. Nakanishi, EPL (Europhysics Letters) 
88 (6), 64001 (2009) 

3. S. Ulrich, T. Aspelmeier, K. Roeller, A. Fingerle, 
S. Herminghaus, A. Zippelius, Physical Review Letters 
102(14), 148002 (2009) 

4. S. Strauch, S. Herminghaus, Soft Matter 8(32), 8271 

( 2012 ) 

5. N. Huang, G. Ovarlez, F. Bertrand, S. Rodts, P. Coussot, 
D. Bonn, Physical Review Letters 94(2), 028301 (2005) 

6 . N. Huang, D. Bonn, Journal of Fluid Mechanics 590, 497 
(2007) 


7. C. Bonnoit, J. Lanuza, A. Lindner, E. Clement, Physical 
Review Letters 105(10), 108302 (2010) 

8 . A. Thornton, D. Krijgsman, R. Fransen, S. Gonzalez, 
D. Tunuguntla, A. te Voortwis, S. Lading, O. Bokhove, 
T. Weinhart, EnginSoft’s SBE&S Magazine 10(1), 48 
(2013). URL http://doc.utwente.nl/89272/ 

9. A. Thornton, T. Weinhart, S. Lading, O. Bokhove, In¬ 
ternational Journal of Modern Physics C 23(08) ( 2012 ) 

10. P. Schall, M. van Hecke, Ann. Rev. Fluid Mech. 42 (2010) 

11. D. Howell, R. Behringer, C. Veje, Physical Review Letters 
82(26), 5241 (1999) 

12. R. Schwarze, A. Gladkyy, F. Uhlig, S. Lading, Granular 
Matter 15(4), 455 (2013) 

13. A. Gladkyy, R. Schwarze, Granular Matter 16(6), 911 
(2014) 

14. S. Herminghaus, Advances in Physics 54(3), 221 (2005) 

15. V. Richefeu, M.S. El Youssoufi, F. Radja’i, in Theoretical 
and Numerical Unsaturated Soil Mechanics (Springer, 
2007), pp. 83-91 

16. P. Darabi, T. Li, K. Pougatch, M. Salcudean, D. Grecov, 
Chemical Engineering Science 65(15), 4472 (2010) 

17. R. Mani, D. Kadau, D. Or, H.J. Herrmann, Physical Re¬ 
view Letters 109(24), 248001 (2012) 

18. R. Mani, D. Kadau, H.J. Herrmann, Granular Matter 
15(4), 447 (2013) 

19. C.D. Willett, M.J. Adams, S.A. Johnson, J.P. Seville, 
Langmuir 16(24), 9396 (2000) 

20. A.M. Alencar, E. Wolfe, S.V. Buldyrev, Physical Review 
E 74(2), 026311 (2006) 

21 . V.A. Lubarda, Acta Mechanica pp. 1-15 (2014) 

22 . F. Soulie, F. Cherblanc, M.S. El Youssoufi, C. Saix, Inter¬ 
national Journal for Numerical and Analytical Methods 
in Geomechanics 30(3), 213 (2006) 

23. A. Singh, V. Magnanimo, S. Lading, Powder Technology 
(2015 submitted) 

24. S. Schollmann, Physical Review E 59(1), 889 (1999) 

25. X. Wang, H. Zhu, A. Yu, Granular Matter 14(3), 411 

( 2012 ) 

26. E. Woldhuis, B.P. Tighe, W. van Saarloos, The European 
Physical Journal E: Soft Matter and Biological Physics 
28(1), 73 (2009) 

27. S. Roy, S. Lading, T. Weinhart, Procedia Engineering 
102, 1531 (2015) 

28. A. Singh, V. Magnanimo, K. Saitoh, S. Lading, New 
Journal of Physics 17(4) (2015) 

29. M. Scheel, R. Seemann, M. Brinkmann, M. Di Michiel, 
A. Sheppard, B. Breidenbach, S. Herminghaus, Nature 
Materials 7(3), 189 (2008) 

30. N. Mitarai, F. Nori, Advances in Physics 55(1-2), 1 
(2006) 

31. A. Denoth, Journal of Glaciology 28(99), 357 (1982) 

32. T. Weigert, S. Ripperger, Particle & Particle Systems 
Characterization 16(5), 238 (1999) 

33. W. Pietsch, R. H, Chemie Ingenieur Technik 39, 885 
(1967) 

34. H. Schubert, Kapillaritdt in pordsen Feststoffsystemen 
(Springer, 1982) 

35. M.A. Erie, D. Dyson, N.R. Morrow, AIChE Journal 
17(1), 115 (1971) 

36. B. Deryaguin. Untersuchungen fiber die reibung und 
adhasion, iv, theorie des anhaftens kleiner teilchen (1934) 

37. G. Lian, C. Thornton, M.J. Adams, Journal of Colloid 
and Interface Science 161(1), 138 (1993) 

38. S. Lading, Particulate Science and Technology 26(1), 33 
(2008) 

39. S. Lading, Particuology 6 ( 6 ), 501 (2008) 



Micro-Macro Transition and Simplified Contact Models for Wet Granular Materials 


15 


40. S. Lading, F. Alonso-Marroqum, Granular Matter 13(2), 
109 (2011) 

41. A. Bouwman, M. Henstra, J. Hegge, Z. Zhang, A. Ingram, 
J. Seville, H. Frijlink, Pharmaceutical Research 22(2), 
270 (2005) 

42. A. Obraniak, T. Gluba, Physicochemical Problems of 
Mineral Processing 48(1), 121 (2012) 

43. G.H. Wortel, J.A. Dijksman, M. van Hecke, Physical Re¬ 
view E 89(1), 012202 (2014) 



